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ABSTRACT 

We examine stochastic variability in the dynamics of X-ray emission from the 
black hole system GRS 1915+105, a strongly variable microquasar commonly 
used for studying relativistic jets and the physics of black hole accretion. The 
analysis of sample observations for 13 different states in both soft (low) and 
hard (high) energy bands is performed by flicker- noise spectroscopy (FNS), a 
phenomenological time series analysis method operating on structure functions 
and power spectrum estimates. We find the values of FNS parameters, including 
the Hurst exponent, flicker-noise parameter, and characteristic time scales, for 
each observation based on multiple 2,500-second continuous data segments. We 
identify four modes of stochastic variability driven by dissipative processes that 
may be related to viscosity fluctuations in the accretion disk around the black 
hole: random (RN), power-law (IF), one-scale (IS), and two-scale (2S). The 
variability modes are generally the same in soft and hard energy bands of the 
same observation. We discuss the potential for future FNS studies of accreting 
black holes. 
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ual (GRS 1915+105) 



INTRODUCTION 



A microquasar is a binary system in which a normal star orbits a compact object, either 
a neutron star or a black hole with a typical mass of 10 solar masses. Under certain condi- 
tions, matter may be transferred from the normal star onto the compact object in the form of 
an accretion disk (for theory, see lShakura fc Sunyaevlll973l ). The subse quent release of grav- 
itatio nal potential energy can lead to the formation of relativistic jets ( iBlandford fc Znajek 
19771 ). as well as strong, vari able emission across a broad range of wavelengths, from ra- 
dio waves to X-ray and 7-ray ( iTananbaum et al.l Il972l ; iMirabel fc Rodriguezl Il994l ) . Due to 
their short characteristic time scales compared to supermassive black holes (i.e., black holes 
having a typical mass of several million solar masses), microquasars are regarded as key 
astrophysical objects for studying relativistic jets and the physics of the accretion process 
f lSams et al.lll996l : IMirabel fc Rodriguezl fl999l McHardv et ~al]l2006h . 



The diverse physical processes in microquasars may be studied indirectly by examining 
temporal variability in emission across the electromagnetic spectrum. Detailed analysis of 
temporal variations may be used to estimate the geometry a nd structure of a binary system 
and to test various astroph ysical models for m icroquasars ( iBellonil |2010| ). This variability 



is commonly characterized (Ivan der Klisll2006l. and references therein) with estimates of the 
powe r spectrum of the ligh t curve (IBelloni et al.ll2002l . 120061 ; iQu et al.ll2010t iMaccarone et al. 



2011 



Neilsen et al.l 120111 ). frequency- dependent time lags between different energy bands 



( Reig et al. 200oh . or a combination of total count r ate and color- color diagrams (colors are 



ratios of X-ray count rates in different energy bands; I B ellonil |2 1 Oh . In the first method, any 
apparent noise components and quasi-periodic oscillations (QPOs) in the power spect ra are 
treated as superpositions of broken power laws and/or Lorentzians ( IBelloni et al.ll2002l ). The 
second method is based on the analysis of phase lag spectra for "soft" (lower) and "hard" 
(higher) energy bands, which are believed to be related to two different sources of X-ray 
emission in the accretion flow; lags between energy bands can help constrain the emission 
geometry. The third method uses ratios between X-ray lightcurves in hard and soft energy 
bands as a proxy for studying the variability of the energy spectrum over time. 

One of the essential features observed in many X-ray emitting astrophysical systems is 
strong apparent randomness on a range of time scales. While the methods above have proven 
successful in discerning and quantifying different variability modes, they could not elucidate 
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the origin of the apparent randomness, which remains an open question. For instance, it 
may be related to variation of some external paramet er, such as stochastic fluctuations of the 
viscosity or mass accretion rate in the accret ion disk (]Lyubarskiilll997l ; iTitarchuk et al.ll2007l ; 
Wilkinson fc Uttley 2009 ; Uttley et alJboil ) ; some complex variabi lity may be due to chaoti c 
(deterministic) fluctuations produced by inner-disk instabilities ( IMisra et al.l I2004L 120061 ). 
Other interpretations are possible, and a complete picture of energy-dependent variability has 
yet to emerge. The origin of apparent randomness in each particular case can be studied using 
the methods of complexity science, such as deterministic chaos, multifractal theory, flicker- 
noise spectroscopy, fractional calculus, cellular automata theory, and other approaches. In 
this paper, we apply a phenomenological complexity science method to analyze one of the 
most remarkable astrophysical systems, the Galactic microquasar GRS 1915+105, which 
consists of a black hole with mass 14±4M , where is the mass of th e Sun, in a 33.5-day 



orbit around a donor star wit 
famous for its re lativistic jets 



h mass 1.2 ± 0.2M a 



different states flBelloni et al. 



Mirabel fc Rodriguez 



2000 



Greiner et al.ll200lf ). GRS 1915+105 is 



1991 and strong X-ray variabili ty in 14 



Klein- Wolt et al. 



20021 ; lHannikainen et all 120051 ). which 



are characterized by highly-structured low-frequency variability in the X-ray lightcurve on 
time scales ranging from seconds to hours, with amplitudes often well in excess of 100% 
(IBelloni et al.l l2000l ). Accordingly, it provides a useful case for studying the underlying 
processes behind the X-ray variability. 

Apparent random behavior in the X-ray li ghtcurves was examined for 12 of the states 



using nonlinear time series analysis techniques (IMisra et al.l I2004L 120061 ; lHarikrishnan et al. 



201ll ). It was suggested that three types of underlying mechanisms may take place: chaotic, 
stochastic (random), and "nonstochastic" (chaotic + colored noise), which implies that sev- 
eral of the 12 variability states could be explained in terms of a deterministic syste m of 
ordinary nonlinear differential equations. The analysis presented by IMisra et al.l ( 120041 ) was 
based on the estimation of correlation dimension, which is complicated in this case by the 
presence of high-frequency Poisson noise in all light curves. Moreover, because the corre- 
lation dimension can approach a constant value for both chaotic data and colored noise, it 
is possible that color ed noise in the light curve da ta for GRS 1915+105 may be incorrectly 
classified as chaotic ( lOsborne fc Provenzald Il989l ). A more advanced procedure including 
surrogate data analysis, singular value dec omposition (SVD) technique, and correlation di- 
mension was applied by IMisra et al.l (120061 ) in an attempt to obtain more conclusive results. 
The analysis of SVD plots for real and surrogate data revealed some apparent attractors 
in the real X-ray variability of GRS 1915+105. However, contamination by noise made it 
difficult to determine whether the attractors were chaotic or periodic (i.e., limit cycles). Cor- 
relation entropy and multi fractal analysis were adde d to the nonlinear time series analysis 
of GRS 1915+105 data by lHarikrishnan et al.l (1201 ll ). This allowed the authors to identify 
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colored noise and estimate the range of fractal dimensions for each observation. It should be 
noted that in addition to the challenges mentioned above, the light curve data are generally 
nonstationary and may exhibit a range of characteristic time scales, which further compli- 
cates the application of nonlin ear time series analysis to these data and makes the results 
inconclusive, as pointed out by lHarikrishnan et al.l ( 120111 ). 

Some of these problems can be overcome by utilizing alternative methods of signal anal- 
ysis for complex systems that are not based on chaos theory. In this study, we apply one such 
approach, flicker-noise spectr oscopy (FNS) ( Timashev et aliboiOa . 2012 ; Timashev fc Polyakov 



20071 ; iTimashevi 12006k 120071 ). a phenomenological method operating on Kolmogorov tran- 



sient structure functions and autocorrelation-based estimates of power spectrum, to ex- 
amine the apparent randomness in the X-ray lightcurves of GRS 1915+105, and compare 
our findings with previo usly reported results obtained by nonlinear time series analysis 
( lHarikrishnan et al.l 1201 if ). In contrast to chaos theory, the FNS method separates between 
regular and stochastic components, handles multiple characteristic time scales, and provides 
a mechanism for excluding the effect of Poisson noise at highest frequencies, thus addressing 
most of the challenges mentioned above. 



2. FLICKER-NOISE SPECTROSCOPY 



Here, we will only deal with the basic FNS relations needed to un derstand the pa- 



rameterization procedure. FNS is described in more detail elsewhere ( Timashev! 



Timashev & Polvakov 


2010a 


)• 



2006 



It should be noted that the term "stochastic" , here and further in this paper, refers to 
random variability in the signals of complex sys t ems characterized by nonlin ear interactions, 
dissipation, and inertia (jTimashev et al.ll2010at Timashev fc Polyakovll2007f ). Conceptually, 
the method separates the analyzed signal into three components: low-frequency regular 
component corresponding to system-specific "resonances" and their interferential (nonlin- 
ear) contributions, stochastic "jump" (random walk) component at larger frequencies cor- 
responding to a dissipative process of anomalous diffusion, and stochastic highest-frequency 
"spike" component corresponding to inertial (non-dissipative) effects. The idea of separat- 
ing the stochastic dynamics of a complex nonlinear system into dissipative random-walk and 
non-dissipative inertial components stems from numerous simulations performed by various 
cellular automata models demonstrating that the stochastic dynamics of these processes is 
associated with intermittency, consecutive alternation of rapid (inertial) changes in the val- 
ues of dynamic variables on small time intervals with small (dissipative) variations of the 
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values on longer time intervals (iBakl 119971 ; IWolframl |2002| ; ISchusterl Il995l ). These studies 
showed that the occurrence of dissipative and inertial changes in such dynamics can be re- 
sponsible for the power law distributions observed in nature, such as Guttenberg-Richter 
and Zipf-Parreto laws and flicker noise. The FNS method further assumes that correlations 
present in sequences of different irregularities, such as spikes and jumps, may be treated 
as main information carriers in the stochastic variability of a complex dissipative system 
( iTimashev fc Polyakovl 120071 ). This assumption has much in common with how wavelet 
analysis evaluates the inverse transform only for those time intervals where the first deriva- 
tive of the signal experiences st eepest changes, i.e., intervals with the highest degree of 
irregular behavior ( jWalnutl 120021 ). The spike-like irregularities are naturally associated in 
FNS with the inertial component, and jump-like irregularities with a dissipative process of 



anoma lous diffusion (random walks), which is illustrated in Fig. 2 of iTimashev fc Polyakov 
( 120071 ) . It should be pointed out that in practical signal analysis, Kolmogorov transient 



structure functions and power spectrum estimates contain only integral contributions from 
stochastic irregularities corresponding to a specific evolution hierarchy level and distributed 
with certain probability densities. If one analyzes the signal at a much higher sampling 
frequency, each of the irregularities observed at a smaller sampling frequency becomes a 
complex structure with its own spike- and jump-like irregularities. Further differences and 
similarities between FNS and other methods for the analysis of the dynamics of complex 
systems, such as nonlinear dynamics (dete rministic chaos), wavelet ana lysis, and fractional 
calculus, are discussed in depth elsewhere (ITimashev fc Polyakovl 120071 ). 



As FNS is a phenomenological method, its main goal is to estimate from analyzed series 
certain model-independent stochastic parameters corresponding to general dissipative and 
inertial processes. These parameters can then be used to develop specific models or refine ex- 
isting ones. The validity of FNS assumptions for many applications was proven by effectively 



geophysical, medical, and p 


rysical systems (" 


Descherevs 


cv et al. 20031; Telesca et al. 


2004; 


Havakawa & Timashev 


2006 




rimashev 2006; 


Ida et al. 


20071 Timas 


rev & Polvakov 


2007; 


Timashev et al. 2009. 5 


>010a 


b 


Rvabinin et al. 


2011: Mirsaidov et al. 


2011; Timashev et al. 


2012h. For instance, it was recently demonstrated bv FNS that the dynamics of X-rav emis- 



sion from GRS 1915+105 and Cygnus X-l, recorded in the timeframe from January 1, 1996 
to December 31, 2005 with a sampli ng; interval of approxim ately 1.5 hrs, manifested a one- 
scale process of anomalous diffusion ( ITimashev et al.ll2010al ). Therefore, it would be logical 
to apply the FNS method to the analysis of the GRS 1915+105 X-ray emission data at much 
smaller (subsecond) sampling intervals, which is the topic of this study. 



In FNS, all introduced parameters for signal V(t), where t is time, are related to the 
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autocorrelation function 

1>(T) = <y{t)v(t + T)) T _ T , (i) 

where r is the time lag parameter (0 < r < T M ), T M is the upper bound for r (T M < T/2), 
and T is the averaging window. The angular brackets in relation (CEJ) stand for the averaging 
over time interval T — r: 

{{■■■)) T -r = j^[~\-)dt. (2) 

The averaging over interval T — r implies that all the characteristics that can be extracted 
by analyzing functions i/j(t) should be regarded as the average values on this interval. 

To extract the information contained in ip(r) {(V(t)) = is assumed), the following 
transforms, or "projections" , of this function are analyzed: cosine transforms (power spec- 
trum estimates) S(f), where / is the frequency, 

[Tit 

S(f) = 2 / (V(t)V(t + t x )) T _ T cos(2 7 r/t 1 ) dh (3) 
Jo 

and its difference moments (Kolmogorov transient structure functions) of the second order 

$(2)( r ) 

& 2 \r) = ([V(t)-V(t + T)} 2 ) T _ T . (4) 

The information contents of S(f) and $( 2 )(r) are generally different, and the parameters 
for both functions are needed to solve parameterization problems. By considering the inter- 
mittent character of signals under study, interpolation expressions for the stochastic compo- 
nents S 8 (f) and 3y C 7 ") of S(.f) an d & 2 Ht), respectively, were derived using the theory of 



generalized functions by iTimashevi (120061 ). It was shown that the stochastic components of 
structure functions $^(r) are formed only by jump-like (random walk) irregularities, and 
stochastic components of functions S(f), which characterize the "energy side" of the process, 
are formed by spike-like (inertial) and jump-like irregularities. It should be noted that r in 
Equations (PP)- (SJ) is considered as a macroscopic parameter exceeding the sampling period 
by at least one order of magnitude. This constraint is required to derive the expressions and 
separate out contributions of dissipative jump-like and inertial (non-dissipative) spike-like 
components. 

The basic idea in parameterization is to use two interpolation expressions for the stochas- 
tic components. The first expression is used to determine the spectral contribution of the 
stochastic components of signal V(t) and exclude the contribution of the low-frequency reg- 
ular component to the parameters related to jump- and spike-like stochastic irregularities. 
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The second interpolation, which deals with the stochastic component in the structure func- 
tion, is used to determine the parameters characterizing the series of jump-like irregularities, 
which correspond to anomalous diffusion (ITimashev et al.ll2010al ). 

For the case of one characteristic scale in the sequences of jump-like and spike-like 
irregularities, the spectral interpolations f or both types of stochastic irregularities have the 
following functional form ( iTimashevi 120061 ) : 



S s (f) 



S,{0) 



(5) 



1 + (2tt/T )-' 

where S s (0), To, and n and are phenomenological parameters, which have different meanings 
for spike- and jump-like irregularities. In other words, the resultant interpolation is a super- 
position of expressions (J5J) for these two components. When the contributions of spikes and 
jumps to the overall stochastic component are comparable in the highest-frequency range 
or the contribution of spikes is negligible, which is true for the data analyzed in this study, 
Equation (JSJ) can be readily used as the resultant interpolation. In this case, S s (0) is the 
phenomenological parameter corresponding to the initial (plateau) value of power spectrum 
estimate, To is the correlation time for jump- and spike-like irregularities after which the 
self-similarity observed in power spectrum estimate breaks down, and n is the flicker-noise 
parameter characterizing the rate of loss of correlations in the series of high-frequency irregu- 
larities in time intervals To. It should be noted that the spectral interpolations for spike- and 
jump-like irregula rities were deriv ed using the theory of generalized functions based on differ- 
ent assumptions (iTimashevi 120061 ) . and generally correspond to different physical processes. 
The negligible contribution of the spike component in this particular case is attributed to 
the high level of Poisson noise at highest frequencies, as will be illustrated later in the paper. 
It should also be pointed out that mathematically Equation (jSJ) represents a superposition 
of flat spectrum at lower frequencies with power law at higher frequencies (broken power law 
with a smooth transition), and T is rel ated to the reciproc al of the characteristic frequency 
used in astrophysical spectral analysis f lBelloni et al.l 120021 ). 



The structure function interpolation for the case of one characterist ic scale in the s e- 
quence of jump-like irregularities (random walks) has the following form ( iTimashevi 120061 ) : 



$< 2 >(r) « 2a 2 [1 - T -\H) ■ T (H, r/T,)} 2 , (6) 

oo 

where F(s,x) = J exp(— t) ■ t 5 " 1 ^ and T(s) = T(s,0) are the complete and incomplete 

X 

gamma functions, respectively (x > and s > 0); a is the standard deviation of the mea- 
sured dynamic variable with dimension [V]; H is the Hurst exponent, which describes the 
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rate at which the dynamic variable "forgets" its values on the time intervals that are less 
than the correlation time T]_. At H = 0.5, interpolation ([6]) corresponds to Brownian motion 
(Fickian diffusion); at H < 0.5, to subdiffusion; and at H > 0.5, to superdiffusion (Levy 
diffusion or Levy flights) (ITimashev et al.ll2010al ). Subdiffusion and superdiffusion are non- 
l inear stochastic processes that are gener ally modeled by fractional Fokker-Planck equations 
(ITsallis fc Lenzil |2002| ; IChen et al.l 1201 0l ). Expression ([6]) at H = 0.5 was shown to corre- 



spond to the diffusion equation with integrodifferential boundary conditions ( ITimashev et al. 
2010ah . 



For the case of two characteristic scales, the resultant spectral interpolation can be 
written as 

S.i(0) S s2 (0) 



Ss(f) 



(2vr/T 01 ; 



nl 



(2tt/T 02 )- 2: 



(7) 



where indices 1 and 2 in subscripts denote respective scales (IParkhutik fc Timashevll2000l ). 
The structure function for two characteristic scales takes the following piecewise form: 



$( 2 ) | 



2oy 



l - 



$i 2) (r) «2<ti 2 + 2a 2 



T(H 2 



Tl2 



r(H 2 



T < T , 



, r > r , 



(9) 



where indices 1 and 2 in subscripts denote respective scales, and To is the value of r cor- 
responding to the transition to the second scale. Expressions ([319]) were derived based on 
the assumption that the characteristic scales are much apart from each other, which justi- 
fies the use of a linear s uperp osition instead of the more complex expression employed by 
Parkhutik fc Timashevl feoooh . 



In some applications, it is also useful to introduce parameters SgiT^ 1 ) (for one-scale 
case) and S^Tq^ 1 ), S^Tq^ 1 ) (for two-scale case) to account f or the "intensity" of ju mp- and 
spike-like irregularities on characteristic frequency intervals (ITimashev et al.ll2012l ). 



FNS interpolations were derived for a wide-sense stationary signal in which the phe- 
nomenological parameters are the same at each level of the system evolution hierarchy. At 
the same time, they can also be applied to the analysis of real signals, which are generally 
nonstationary, but can be characterized by a standard deviation within a specific averaging 
window. In this case, the real signals at specific averaging intervals and sampling frequen- 
cies should be regarded as quasi- stationary with certain values of standard deviation and 
other phenomenological parameters. It should be noted that the values of phenomenological 
parameters may vary on different quasi-stationary intervals. 
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The detailed FNS parameterization algorithms for one- and two-scale cases in discrete 
form, which were adapted for the X-ray emission data, are described in Appendices A and 
B, respectively. 



3. RESULTS AND DISCUSSION 

In order to sample the stochastic variability in GRS 1915+105, we analyze 15 of the 
observations of this system (13 different X-ray states) made with the Rossi X-ray Timing 
Explorer over the last 15 years, which are illustrated in Figure [TJ Each observation consists 
of several ~ 2, 500-second intervals separated by gaps. The specific observations are selected 
to study at least one long instance of each of GRS 1915+105's known states (except for the cu 
state, for which we did not find a long enough continuous data set). From each observation, 
we extract X-ray lightcurves in two energy bands (soft: 2-6 keV, and hard: 15-60 keV) with 
sampling frequencies fa of 16 and 128 Hz. We also extract X-ray background lightcurves 
in both energy bands to subtract from each time bin, and we correct the lightcurve count 
rates for detector dead time. Taking into consideration that all analyzed time series exhibit 
Poisson noise in the highest frequency range, we add a term to our interpolations (IMH]) 
accounting for this noise. For the relatively low frequencies considered here, the expre ssion 



for the dead-time-corrected Poisson noise level (jZhang et al.lll995l ; iMorgan et al.l 119971 ) can 
be approximated by a constant value (in general, it is frequency- dependent; see Appendix 
C). Considering only continuous 2,500-s intervals in two energy bands, we analyze 94 time 
series in total. 

We perform the analysis at Tm = T/9 ~ 280 s to make sure the autocorrelation-based 
power spectrum estimate is close to the actual power spectrum value. The spectral FNS 
parameters a re also found for the power spectrum estimate calculated by the Welch method 



(I Welch! 1 196 71 ) (averaging 8 windows of size T/9 with a 50% overlap). The Welch method 



i s similar to metho ds for computing power spectra commonly used in X-ray astrophysics 



( IMorgan et al.lll997l ). The difference between the parameter values found by the two methods 
rarely exceeded 10%, which implies that the autocorrelation-based power spectrum estimate 
is adequate for our analysis. 

For the 15 observations analyzed in this study, the variations of FNS parameters within 
a given observation were not significant in most cases (generally less than 10%), which 
suggests that the parameterization method is reliable and that specific variability modes 
last longer than the observations studied here. Table 1 lists the observation-averaged values 
of parameters for the 12 observations where the contribution of Poisson noise is significant 
in the frequency range close to 16Hz. The other 3 observations, which have a relatively 
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Fig. 1. — First 2,500-second record of continuous activity for light curves of 15 sample 
observations (soft X-ray energy band). 
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high signal-to-Poisson-noise ratio in that frequency range, are parameterized for a higher 
sampling frequency fd = 128 Hz (Table 2). The parameters S^Tq" 1 ), S^T^ 1 ), S a (T^), a, 
ax, and cr 2 are not listed in the tables because their values cannot be reliably determined 
when a significant Poisson noise level is present. 



Table 1: Parameters for the light curves of 12 GRS 1915+105 observations with IS, IF, or RN behavior (/^ = 16 Hz) 









Soft X 


-rays 








Hard X-rays 


Observation ID 


Class 




Hv 2 , 1/s 


Type 


PN% 3 


n 


H 


To, s 


Ti, s 


Hv 2 , 1/s 


Type 


PN% 3 


n 


H 


To, s 


Ti, s 


90105-06-01-00 


7 


4 


2790 


IS 


15 


1.73 


0.37 


13 


17 


122 


IS 


30 


2.15 


0.59 


5 


6 


40703-01-24-00 


K 


3 


1443 


IS 


20 


2.55 


0.77 


■5 


5 


89 


IS 


35 


2.61 


0.86 


3 


3 


60405-01-02-00 


P 


3 


1548 


IS 


15 


2.15 


0.58 


6 


7 


100 


IS 


35 


1.91 


0.47 


2 


3 


20402-01-54-00 


6 


3 


2134 


IF 


15 


1.2 


0.1 






160 


IF 


■50 


1.04 


0.06 


- 


- 


20402-01-37-01 


A 


2 


1409 


IF 


25 


1.91 


0.46 






182 


IF 


30 


1.57 


0.29 






92092-03-01-00 


/' 


4 


3814 


IF 


20 


1.7 


0.35 






148 


IF 


45 


1.94 


0.48 






80127-01-02-00 




4 


889 


IF 


10 


1.47 


0.25 






44 


IF 


65 


3.06 


0.92 






30183-02-02-02 




3 


611 


IF 


40 


1.36 


0.21 






49 


IF 


70 


1.88 


0.3 






20187-02-01-00 


a 


2 


581 


IS 


10 


1.53 


0.27 


50 


86 


129 


RN 












50125-01-04-00 


V 


3 


1102 


IF 


5 


1.29 


0.16 






107 


RN 












50405-01-01-00 


V 


2 


848 


IF 


5 


1.37 


0.18 






92 


RN 












20402-01-13-00 


X 


3 


387 


RN 












119 


RN 













to 



1 number of intervals in observation 

2 average count rate 

3 percentage of Poisson noise in the highest-frequency range - determined for the plot in double logarithmic 
scale 

Table 2: Parameters for the light curves of 3 GRS 1915+105 observations with 2S behavior (/^ = 128 Hz; notes same 
as in Table 1) 



Observation Id 


Class 




Band 


(iv, 1/s 


Type 


m 




T i,s 


T n ,s 


712 


H 2 


T]2,S 


2l2,8 


92092-01-01-00 




3 


Soft 


1604 


2S 


1.73 


0.53 


0.22 


0.20 


3.07 


1.06 


12 


10 








Hard 


121 


2S 


1.76 


0.58 


0.16 


0.14 


2.49 


0.66 


30 


36 


20186-03-02-01 





4 


Soft 


1738 


2S 


2.53 


0.88 


0.13 


0.13 


2.3 


0.65 


34 


42 








Hard 


228 


2S 


2.19 


0.75 


0.11 


0.10 


2.5 


0.7 


38 


46 


20402-01-45-02 


e 


4 


Soft 


1783 


2S 


2.38 


0.82 


0.13 


0.13 


2.21 


0.58 


39 


56 








Hard 


232 


2S 


1.9 


0.66 


0.13 


0.11 


2.68 


0.74 


35 


44 



13 



Our analysis shows that the noise behavior for the observations considered in this study 
exhibits four different modes of stochastic variability: "random" (RN), power-law (IF), 
one-scale (IS), and two-scale (2S). The RN variability, which is illustrated in Fig [2(a), 
corresponds to values of n less than unity and the overall variation of spectral power in 
the considered frequency range by less than two orders of magnitude. The IF variability is 
characterized by a power-law (l// n ) dependence at low and intermediate frequencies and a 
flat Poisson noise level at highest frequencies, as can be seen in Figure [2(b). This implies 
that the stochastic process for this variability mode exhibits a form of self-similar behavior 
and is far from reaching the steady state (characteristic time scale) at the averaging interval 
and frequency range considered here. The IS variability follows the interpolation (jSJ) at 
low and intermediate frequency ranges and exhibits a Poisson noise at highest frequencies, 
as illustrated in Figure [2(c). The 2S variability is adequately described by interpolation 
flZJ) and contains two characteristic time scales, as can be seen in Figure [2(d). In Figure 
3, we present the structure functions (panels (a) and (c)) and their stochastic components 
(panels (b) and (d)) for sample IS and 2S variability modes, respectively. The stochastic 
component is obtained by subtracting the resonant component from the overall structure 
function. For the IS mode, the stochastic com ponent is described by an anomalous diffusion 



process reaching a steady state (Figure [3(b); iTimashev et al.ll2010al ). while for the 2S case 
a two-scale process is implied (Figure [3(d)). 

It can be seen that the parameters for two instances of the same class {y in Table 
1 and 9 in Table 2) stay approximately same, which suggests that a particular stochastic 
variability (certain variability mode and ranges for parameter values) may characterize a 
class of observations. However, an exhaustive analysis of all X-ray emission observations for 
GRS 1915+105 is needed to draw definite conclusions in this context. 

For most of analyzed intervals corresponding to variability modes IF, IS, and 2S, we 
noticed t hat n ~ 2H + 1 (n-\ ~ 2H^ + 1, n 2 ~ 2H 2 + 1), as expected for fractional Brownian 



motions (iMalamud fc Turcottdll999l ). We also observed that the characteristic time scales 
in power spectrum estimate and structure function were close to each other for most of 
the intervals characterized by variability modes IS and 2S. Both of these facts imply that 
the random-walk component plays the dominant role in the stochastic variability of X-ray 
emission, and the contribution of highest-frequ ency inertial spike-irregularities in these data 



is not significant (ITimashev fc Polyakovl 120071 1 . This can most likely be attributed to the 



"erasing" effect of Poisson noise at highest frequencies, which is illustrated in Figure [2] and 
by the values of parameter PN% in Table 1. 



Tables 1 and 2 as well as interval-specific analysis suggest that there is no statistically 
significant difference between the FNS parameter values for soft and hard X-ray energy bands 
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f[Hz] f[Hz] 

Fig. 2. — Typical power spectrum plots for each mode of stochastic variability [thick black 
line is spectral FNS interpolation; thin gray line is experimental power spectrum estimate]: 

(a) RN, class x> observation 20402-01-13-00, l s ^ interval, soft X-ray energy band, n m 0.52; 

(b) IF, class fi, observation 920920-03-01-00, 3 rd interval, soft X -ray energy band; (c) IS, 
class 7, observation 90105-06-01-00, 3 rd interval, soft X-ray energy band; (d) 2S, class 9, 

th 

observation 20402-01-45-02, 4 lli interval, soft X -ray energy band. 
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Fig. 3. — Sample structure function plots for variability modes IS and 2S [1 (blue)- experi- 
mental curve, 2 (green) - interpolation, 3 (red) - resonant component of structure function]: 
(a) structure function plot for IS, class 7, observation 90105-06-01-00, 3 rc ^ interval, soft X-ray 
energy band and (b) its stochastic component of structure function; (c) structure function 
plot for 2S, class 9, observation 20402-01-45-02, 4 U1 interval, soft X -ray energy band and (d) 
its stochastic component of structure function. See the electronic edition of the Journal for 
a color version of this figure. 
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in the same interval. The only differences are in the count rate magnitude (it is smaller for 
hard X-rays) and the frequency range dominated by the Poisson noise (PN% is higher for 
hard X-rays). Soft X-ray emission is nominally associated with the thermal spectrum of 
the accretion disk around the black hole, while hard X-ray emission is generally believed 
to be related to Compton scattering in a hot electron corona around the disk. The fact 
that the FNS parameters for these energy bands are close in value could suggest that at 
these frequencies, there is a common mechanism producing stochastic variability. Future 
analysis in more energy bands could help determine whether the FNS parameters are truly 
energy- independent . 

Let us compare the classification in Tables 1 and 2 with the results obtained by the 
nonlinear time seri es analysis of total count rates for sample observations corresponding to 
12 different states (IHarikrishnan et al.ll201ll . Table 1) . As total count rat e s are dominated by 
the soft X-ray energy band, we compare the results of lHarikrishnan et al.l (120111 ) with the soft 
X-ray classification presented in this paper. The noise variabilit y for class y is chara c terize d 
as random by both analyses. Class <fi is marked as random by IHarikrishnan et al.l ( 120111 ). 
but is assigned to IF here. The difference, in our opinion, is due to the significant effect 
of the Poisson noise within half of the considered frequency range (high values of PN%), 
which was not separately taken into account i n the nonlinear time series analysis. Class 7 
was marked as random by IHarikrishnan et al.l (120 111 ), but is categorized as IS by FNS. This 
discrepancy could be d ue to the plateau in the l ow-frequency range of the power spectrum 
(Figure [2] in our study; IHarikrishnan et al.l 1201 ll . Figure 10): although it may mimic white 
noise in the nonlinear time series analysis, in FNS the plateau corresponds to the loss of 
correlations ( at r > T n ) within a s equen ce of stochastic irregularities. Class 5 is categorized 
as random by IHarikrishnan et al.l ( 120111 ). but is marked as IF in FNS. Considering that the 
value of n is close to unity in this case (boundary range), additional sample observations 
for this class would be needed to determine which method is more accurate. Classes k, 



A, and fi were labeled by IHarikrishnan et al.l (120111 ) as "deterministic nonlinear + colored 
noise". All three are marked as IF or IS in this study, which corresponds to a variant of 
"colored noise". In classes where there are strong periodic or quasi- regular ("resonant") 
components (e.g. a, (3, 6, v, k, A, and p \ it is difficult to compare our categorizations 
directly to the results of iMisra et al.l (120061 ); IHarikrishnan et al.l (120111 ). because in contrast 
to those studies, we subtract the resonant components prior to performin g our analysis. In 



other words, the regular (deterministic) behavior observed in SVD plots ( IMisra et al.l 12006 



Harikrishnan et al.l 1201 11 ) could be dominated by the resonant periodicities, rather than the 
stochastic component studied here. Moreover, classes j3 and 6 contain two separate stochastic 
scales in different frequency ranges, which would be difficult to capture in terms of chaos 
theory, originally developed for systems of nonlinear ordinary differential equations with one 



-17- 



characteristic scale. 



The above analysis demonstrates that the random-walk component plays the dominant 
role in the stochastic variability of modes IF, IS, and 2S. The random- walk component 
can be interpreted in terms of a dissipative process of anomalous diffusion ( iTimashev et al. 
2010al ). It has been argued that stochastic fluctuations of the viscosity in the accretion 
flow may be responsible for X-ray variability on a wide range of t ime scales, as well as 
the presence of multiple characteristic time s cales (jLyubarskiil 119971 ; iTitarchuk et all 120071 ; 
Wilkinson &: Uttleyl l2009t lUttley et all 1201 ll ). This is because variations in the viscosity 
lead to changes in the accretion rate, which propagate towards the black hole, resulting in 
flickering X-ray emission. In principle, the radial dependence of the amplitude of the viscous 
fluctuations can determine the slope of the power spectrum, and different components in the 
accretion flow (e.g. corona, inner accretion disk, and outer accretion disk) can contribute at 
different frequencies. The stochastic modes described here may be related to these viscous 
fluctuations in the accretion disk, in which case a more comprehensive analysis with FNS 
could provide a useful characterization of X-ray states from the perspective of stochastic 
variability. 



4. CONCLUDING REMARKS 

Our analysis shows that the consideration of different frequency ranges in flicker-noise 
spectroscopy allows one to analyze the stochastic variability in complex signals of GRS 
1915+105, despite the significant contribution of Poisson noise in the highest frequency 
range, high low-frequency quasi-regular variability, and presence of multiple stochastic time 
scales. As a proof of concept for future study, the examination of sample observations for 
13 different states demonstrates that there are at least four types of stochastic variability in 
the X-ray emission of GRS 1915+105: random (RN), power-scale (IF), one-scale (IS), and 
two-scale (2S). The last three are related to random walk processes interpretable in terms of 
a dissipative process of anomalous diffusion. In this case, the random walk processes could 
be related to viscosity fluctuations in the accretion disk. 

Other stochastic variability modes are also possible in the frequency ranges under con- 
sideration. For example, it is natural to expect a stochastic variability mode representing 
an intermediate category between IS and 2S, i.e., a superposition of one full scale (IS) and 
self-similar (IF) behavior in different frequency ranges. This is borne out by astrophysi- 
cal observations, which exhibit a wide range of extremely complex po wer spectra, and are 



frequently modele d as the sum of a number of different components (jBelloni et all 12002 



van der Klisl 120061 ) . An exhaustive analysis of all available X-ray emission observations for 
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GRS 1915+105 is needed to get a complete list of variability modes, obtain more accurate 
values of FNS parameters for each known class of the binary system, and revisit the estab- 
lished classification from the stochastic perspective. A similar analysis of other accreting 
black holes could help to identify the effect of the extreme variability of GRS 1915+105 on 
the FNS parameters, and to determine the true relationship of these parameters to stochas- 
tic processes in the disk. It would also be interesting to study the va riability at other 
wave l engths, e.g., in obs ervations of the radio jet with the Ryle Telescope (jPooley fc Fender 



19971 ; Fender et al.lll997l ) or the EVLA, via FNS parameterization and cross-correlation with 



X-rays, in order to compare stochastic variability in different accretion processes. 
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A. Parameterization algorithm for one-scale case in discrete form 

Consider a time series Vd(k). The subscript d here and below is used to denote the 
discrete form of expressions. Let N t be the number of points corresponding to the selected 
averaging interval T, M be the number of points used in estimating the autocorrelation 
function. In this case, the parameterization procedure can be written as follows: 

1. Calculate the arithmetic mean for the signal: 

1 Nt 

* fc=i 

2. Subtract the arithmetic mean from the series V d (k): 

V d (k) = V d (k)-fi v . (A2) 

3. Calculate the autocorrelation function for the series V d : 

Nt-p _ 

MP) = tt^ E V d (k) V d (k + p),p = 0..M. (A3) 
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The autocorrelation interval M should not be higher than N t /A (higher values of M will 
result in the loss of statistical information in estimating the autocorrelation function). To go 
from discrete form to the continuous one, one can use the following expression: p = N t r/T. 

4. Calculate the discrete cosine transform of the autocorrelation function: 

S d (q) = MO) + fa (M) (-1)" + 2 MP) cos (^f), (A4) 

p=i 

where q = 0..M. For q = 1..M — 1, Sd(q) should be multiplied by 2, which is the standard 
procedure for discrete Fourier transforms to take into account the spectral values in the 
second half of the frequency range. Here, relations q = 2ff d ~ 1 M and Sd{q) = S(f) x 
fd describe the equivalence between the discrete and continuous forms of power spectrum 
estimate. 

5. Calculate S s d(0) as the average value of the power spectrum for the points 1 and 2 (point 
0, which corresponds to the zero frequency, is not used in calculating S s d(0)): 

Ssd (0) = . (A5) 

6. Interpolate |Sd(g)| given by Equation (1A4|) using the expression: 

S sd (q) = ,J S fl , n +Ps (A6) 

by the method of nonlinear least-square fitting to determine the values of parameters n and 
T d- The constant term P$ is an estimate of the Poisson noise level calculated as the average 
value of Sd for 100 highest- frequency points. The fitting is done on the basis of a double 
logarithmic scale, dividing the entire series into a set of equal intervals (in the calculations 
presented in this study we took 200 equal intervals). We used the t rust-region algorith m for 



nonlinear square fitting, which is built in MATLAB v. 7 or higher ( [Branch et al.lll999[ ). 
7. Separate out the resonant component: 

S rd (q) = S d (q) - S sd (q) , q = 0..M. (A7) 



8. Calculate the autocorrelation function for the resonant component as the inverse discrete 
cosine transform of S r d(q). When q = 1..M — 1, divide S r d(q) by 2 to take into account the 
spectral values in the second half of the frequency range. Then calculate the inverse cosine 
transform: 

Vw(p) = 7^7 {Srd(0) + S rd (M) (-1H 




M-l 



2 cos \. (A8) 



M 
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9. Calculate the difference moment for the resonant component: 

$S(P) = 2 [Vrd(0) - <Mp)] ,p = 0..M. 

The continuous equivalent of & r J(p) is $r (t). 

10. Calculate the difference moment for the experimental series: 



1 



P 



N t -p 

E 

fc=i 



V d {k)-V d (k + p) 



(A9) 



(A10) 



(All) 



11. Calculate the difference moment for the random component: 
The continuous equivalent of & esd (p) is $es (t). 

(2) 

12. Determine the parameters a, if, Tid, P$ by fitting $ c J d (p) in Equation (]A11|) to the 



interpolation expression of the anomalous diffusion type (ITimashev et al.ll2010ar ): 



^3(p) = 2a 2 x [1 - T-\H) x T(H,p/T ld )Y + P 9 , 



(A12) 



where T(s,x) = J exp(— t)t s 1 dt, T(s) = r(s, 0), using the same least-square fitting method 

X 

as in step 6. Here, the constant term P$ is an estimate of the Poisson noise contribution to 
the structure function. 

13. Calculate S sd (T^ d ) by Equation ( 1A6[) . 

14. After the values of all six FNS parameters - a, T od , T ld , H, n, S s d(T^) - are determined, 
calculate the dimensional values for T od , T ld , S sd (TQ d ): T = T od x At, T x = T xd x At, 
S s (Tv l ) = SsdiT-J) x At, where At = . 

(2) 

15. Calculate the relative error e$ in the interpolation of difference moment $^ (p): 



M 

E 

P =i 



(p) - ^ (p) - C j (p) 



,(2) 



(2) 



M (?) 

E $i 2) (P) 

p=l 



x 100%. 



(A13) 



Here, the error is determined as the ratio of the difference of areas between the experimental 
structure function and the total interpolation function to the area of the experimental struc- 
ture function. The areas are calculated by numerical integration using the rectangle method 
because the original se ries have a rather large number of points. The parameterization is 
successful if e$ < 10% (ITimashev et al.ll2010al ). 
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B. Parameterization algorithm for two-scale case in discrete form 

Steps 1 through 5 are same as in the one-scale case (Appendix A). 
In step 6, the following interpolation is used 

q S Bld (0) s sd (Q) - s ald (o) , p mn 

Here, the method of nonlinear least-square fitting is applied in log-log scale to determine the 
values of S sld (0), T Qld , m, T 02d , and n 2 . 

Steps 7 through 11 are same as in the one-scale case (Appendix A). 

In step 12, we first use the following interpolation to estimate the values of parameters 
for the first scale: 



$iJ(p) = 2^x[l 



T(H x ,p/T lld ) 



0"! 



r(#i 

T(H 2 ,p/T 12d )\ 2 



r(# 2 



+ p*, 



(B2) 



The nonlinear least-square fitting method in log-log scale is applied to find the values 
of parameters cri, Hi, T nd , a 2 , H 2 , and Ti 2d . The parameter P$ is estimated as the value of 
^esd(p) a ^ P = The se ^ corresponding to the lower value of T ud and T 1M is used as the 
values for o\, Hi, and Tn d . 

The parameter rod is estimated as the value of p after which $2 d (f>) — P$ exceeds 2crf 
by at least 3%. Then we use interpolation 



*2(p) = 2*? + 2a a a 



1 - 



r(# 2 , ^ 



1 2 



T(H 2 



P* 



(B3) 



for p > r od to find the values of a 2 , H 2 , and T\ 2d . 

At step 13, calculate S sd {T^ ld ) and ^(Pqm) by Equation f lBlj) . 

At step 14, calculate the dimensional values for P ow , P 02d , Piw, Pi2d, ^(P^), S sd (T~ 2d ): 
Poi = P ow x At, P 02 = P 02d x At, P n = P 1W x At, P 12 = P 12d x At, ^(P^ 1 ) = SW^) x At, 
Ssi^ 1 ) = S sd (T- d ) x At, where At = f, 1 . 



Step 15 is same as in the one-scale case (Appendix A). 
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C. Approximation for dead-time-corrected Poisson noise level expression 



The formul a for computing; the Poisson power cor rected for instrument dead time is 
given in Refs. fjZhang et al.l Il995t iMorgan et al.lll997l ). Adapting this for power spectra 
based on the autocorrelator, we have: 



S P (f) = A x (2 



2r r d ( 1 - 



r — r oT d ) 



+2r vle r Q rf le 



Sin (TTTyief) 



(CI) 



where A is a scale factor between FNS units and standard astrophysical power spectra (where 
pure Poisson noise has Sp(f) = 2), r is the count rate per proportional counter unit (PCU), 
N t is the number of frequency points, r v \ e is the rate of very large events (e.g. cosmic rays; 
VLE) per PCU, U is the bin size, is the dead time per event, r„; e is the dead time window 
for each VLE. In our case, r vie = 200 counts/s/PCU, tfe=0.0625 or 0.0078 s, r d = 10/xs, and 

T v i e = 150/iS. 

It is more convenient for analysis to write Equation (ICljl as 



S P (f) = B + Ccos [Df] + E, 



(C2) 



where 
B 

C 
D 



2A 



1 - 2r r d 1 



2L) 

2r 6 J 



2irt h 



E = 2Ar vle r r 1 



vie 



sin(7rr„ ;e /) 



2Ar v i e r T? 



vie' 



It can be easily shown for our data that B ^> C + E (usually by at least 2 orders of 
magnitude). This implies that in our study we can use the constant- value approximation 



Sptf) w B. 



(C3) 
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